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ABSTRACT 


~A  graphical  test  bed  in  which  the  results  of  a  simulation 
experiment  can  be  reported  and  analyzed  is  described.  The 
test  bed  is  based  on  the  regression  adjusted  graphics  and 
estimation  (RAGE)  methodology  developed  by  Heidelberger  and 
Lewis  p^f.^-1]  for  regenerative  simulations.  From  the  graphics 
and  the  associated  numerics  the  experimenter  can  summarize 
and  see  simultaneously  relative  properties,  such  as  bias, 
normality  and  standard  deviation,  of  several  estimators  of  a 
characteristic  of  a  population  for  up  to  eight  sample  sizes. 

The  graphics  is  supported  on  a  line  printer  to  make  it  and 
the  program  portable. 
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I.  INTRODUCTION 


During  the  course  o£  a  simulation  experiment  the  experi¬ 
menter  is  frequently  confronted  with  the  problem  of  estimating 
some  function  or  parameter  from  the  generated  data.  He  will 
usually  have  a  choice  of  one  or  more  estimators,  parametric 
or  nonparametric ,  that  can  be  used  to  estimate  the  function 
or  parameter.  These  estimators  will  vary  in  their  statistical 
properties.  Some  will  estimate  with  low  variance  but  contain 
bias,  others  will  have  just  the  opposite  characteristics  and 
still  others  will  contain  varying  degrees  of  both.  How  these 
different  estimators  compare  in  relation  to  each  other  as 
well  as  over  varying  sample  sizes  are  of  concern  to  the  experi¬ 
menter.  Comparison  of  estimators  are  important  in  that  they 
help  the  experimenter  to  choose  the  most  appropriate  estimator 
for  his  particular  analysis.  Comparisons  may  also  be  helpful 
in  deciding  whether  a  new  estimator  is  of  any  use.  The  appro¬ 
priate  estimator  to  use  is  a  function  of  the  bias,  variance 
and  nonnormality  in  the  estimate  for  available  run  length 
(sample  size). 

When  the  simulated  data  from  which  the  properties  of  the 
estimators  is  to  be  examined  is  independent  and  identically 
distributed  (iid)  it  can  be  sectioned  into  M  independent 
blocks  of  specified  sample  size  (N) .  From  each  of  these 
blocks  the  estimator  can  be  applied  to  all  subsamples  of 
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size  rij^,  thereby  creating  samples  of  estimates  of  varying 

size  for  each  estimator.  As  a  first  cut  at  the  properties 

of  the  estimators,  box  plots  of  these  samples  of  estimates 

are  drawn  for  subsample  sizes  n  (see,  e.g.  Figure  3a). 

k 

However,  something  is  usually  known  about  the  properties 

of  the  estimators  and  these  properties  can  be  used  to  enhance 

the  graph  of  box  plots  versus  subsample  size,  nj^.  Thus  it 

is  often  known  that  bias  in  the  estimator  is  a  polynomial 

in  1/n,  and  the  variance  of  the  estimator  is  a  polynomial  in 

1  ,  1  ,  1  .  . . .  Also  the  estimators  are  usually  asymp- 

n  ^  n^ 

totically  normally  distributed.  Consequently  the  graphs  are 
enhanced  by  an  estimated  unbiased  regression  line  for  the 
mean  of  the  estimates  of  the- parameter  as  well  as  a  horizontal 
line  for  the  regression-estimated  expected  value  of  the 
parameter  (asymptote).  A  regression  is  also  performed  for 
the  variance  structure  of  the  estimators  and  the  sample 
characteristics  (kurtosis  and  skewness)  are  given  under  each 
boxplot  to  assist  in  judging  the  convergence  to  normality. 

Heidelberger  and  Lewis  [Ref.  1]  used  this  methodology  in 
the  context  of  regenerative  simulation  to  develop  a  procedure 
and  protocol  to  produce  unbiased  and  normally  distributed 
regression  adjusted  regenerative  estimates  (RAGE) ,  with  a 
given  precision  or  for  total  available  computing  time.  How¬ 
ever,  the  graphics  associated  with  this  procedure  required  a 
Tektronic  4013  terminal  using  an  APL  graphics  package  making 
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its  use  limited  to  those  organizations  possessing  such  a 
terminal.  The  protocol  also  requires  a  relatively  experi¬ 
enced  user  to  properly  interpret  the  graphical  and  statisti¬ 
cal  output  and  to  apply  the  sequential  steps  of  the  protocol 

The  purpose  of  this  thesis  is  to  apply  the  methodology 
developed  by  Heidelberger  and  Lewis  [Ref.  1]  to  comparing 
the  bias  and  variance  of  estimators  of  a  parameter  from  simu 
lated  output.  This  objective  also  includes  supporting  the 
graphics  associated  with  the  methodology  on  a  line  printer 
to  make  its  use  more  accessible. 

Details  of  the  graphics  and  methodology  are  given  in 
Chapter  II,  entitled  USAGE.  Details  of  the  arguments  of  the 
RAGE  program  are  given  in  Chapter  III  and  a  sample  program 
using  the  RAGE  graphics  is  given  in  Chapter  IV.  All  the 
user  has  to  supply  is  the  input  data  and  FORTRAN  functions 
to  compute  the  estimators  for  any  given  subset  of  data. 

A  comparative  study  of  estimators  of  the  gamma  shape 
parameter  using  the  graphical  test  bed  RAGE  is  given  in 
Chapter  V.  This  comparative  study  illustrates  the  different 
properties  of  estimators  that  can  be  seen  with  regression 
adjusted  graphics  as  sample  size  increases. 
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1 1 .  USAGE 


The  regression  adjusted  graphics  and  estimating  (RAGE) 
program  is  designed  to  evaluate  distributional  properties  of 
one  or  more  estimators  as  a  function  of  sample  size  from  a 
stationary  series  of  raw  data.  These  estimators  are  program¬ 
med  as  functions  by  the  user  and  declared  external  in  the 
main  program.  The  graphical  and  numerical  results  that 
accompany  each  estimator  are  varied  and  will  depend  to  some 
extent  on  the  options  desired  by  the  user.  The  purpose  of 
the  following  discussion  is  to  explain  the  graphical  and 
numerical  output  as  well  as  how  they  were  obtained  and  some 
of  the  options  available  to  the  user. 

The  general  procedure  begins  with  a  given  iid  data  base 
or  one  that  has  been  simulated  and  an  unknown  parameter  or 
population  characteristic,  0,  that  is  to  be  estimated  by  an 

A 

estimator,  9,  from  the  data.  The  experimenter  must  first 
decide  on  the  number  of  independent  blocks  or  sections  (M) 
of  size  N  the  data  is  to  be  sectioned  into.  These  M  sections 
or  blocks  will  henceforth  be  referred  to  as  replications; 
the  terminology  will  seem  more  appropriate  as  the  discussion 
continues.  The  number  of  replications  (M)  can  alternatively 
be  determined  indirectly  by  the  sample  size  (N)  desired  in 
each  section.  In  either  case  MxN  must  not  exceed  the  total 
allowable  sample  size.  This  sectioning  of  the  data  into 
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M  replications  of  sample  size  N  is  represented  in  Figure  1 
where  the  entire  block  represents  the  MxN  samples  from  the 
total  sample  population. 


Within  each  replication  the  experimenter  must  decide  on 

A 

the  subsample  sizes  that  the  estimator,  0,  is  to  use  in 
estimating  the  parameter.  The  number  of  different  subsample 
sizes,  K,  that  can  be  used  to  resection  the  same  data  can  vary 
from  one  to  eight  and  the  size  of  each  subsample  can  vary 
from  one  to  the  size  of  the  sample,  N.  This  repeated  sub¬ 
sectioning  of  the  same  data  is  represented  pictorially  by 
Figure  2  where  [x]  refers  to  the  highest  integer  less  than 
or  equal  to  x  and  rj^,  k=l ,  ...,  K  refers  to  the  number  of 
estimates  possible  from  a  sample  of  size  N  using  subsample 
sizes  njj,  k“l,  ...»  K  to  estimate  the  parameter.  Each  box 
or  rectangle  in  Figure  2  represents  the  same  MxN  data  set 
with  the  same  M  replications  of  size  N  subsectioned  into 
different  sections  of  size  nj^.  Note  that  the  nj^'s  divide 
N  evenly  in  Figure  2.  This  is  not  nece-ssary  but  certainly 
desirable  since  all  available  information  will  be  utilized. 


Sample 

size 


Total  Samples  MxN 


Population  of 
Section  Estimates  for 

Rep.)  subsample  size  n 


population  of  estimates  from  each  subsample  •••>  K) 

box  plots  are  plotted  so  that  the  distribution  of  the  esti¬ 
mates  for  the  varying  subsample  sizes  can  be  seen  graphically. 

In  addition  to  these  box  plots  of  the  estimates,  a 
regression  adjusted  estimate  (RARE)  of  the  parameter  is 
computed  and  plotted  over  the  boxplots.  In  addition  to  the 
RARE  estimate,  unbiased  RARE  estimates  out  to  a  specified 
number  of  terras  determined  by  the  user  are  computed  and  plotted 
on  the  graph  (The  mathematical  derivation  of  the  RARE  and 
unbiased  RARE  estimates  are  given  later  in  this  section) . 

The  unbiased  RARE  estimates  are  computed  as  a  function  of 
subsample  size  and  plotted  in  conjunction  with  the  RARE 
estimate  so  that  a  comparison  of  the  two  values  can  be  made 
graphically  and  the  evolution  to  an  unbiased  state  can  be 
seen  as  the  subsample  size  increases. 

This  idea  can  be  made  clearer  through  the  use  of  an  example 
and  its  associated  graphics  and  numerics.  Figure  4a.  Here  an 
independent  and  identically  distributed  sample  of  50,000 
gamma  distributed  random  numbers  (shape  parameter,  K=S.O) 
has  been  generated  and  sectioned  into  M*100  sections  (replica¬ 
tions)  with  each  replication  containing  N»500  randon  numbers 
(sample  size).  This  information  is  printed  at  the  top  of 
the  graph  for  future  identification.  The  experimenter  chose 
eight  different  subsample  sizes  at  which  to  examine  the 
parameter  (33,  50,  71,  125,  166,  250  and  500).  These 
subsample  sizes  (n^^'s)  are  printed  on  the  horizontal  scale 


of  the  graph.  The  population  of  estimates  of  the  parameter 
at  each  subsample  size  are  then  plotted  as  vertical  box  plots 
above  these  numbers.  When  the  population  of  estimates  deter¬ 
mined  from  a  subsample  size  is  less  than  nine,  box  plots 
are  omitted  and  the  values  of  these  estimates  are  simply 
plotted  with  the  sumbol  "0"  on  the  graph.  There  are  no 
estimates  with  populations  totaling  less  than  nine  in  Figure 
4a  but  the  graphical  output  for  such  cases  can  be  seen  in 
Figure  3a  where  subsaraple  sizes,  125,  166,  250  and  500  total 
less  than  nine  estimates. 

The  horizontal  scale  of  the  graph  is  scaled  roughly  to 
the  minimum  and  maximum  subsample  sizes  with  the  box  plot  of 
each  subsample's  estimates  spaced  proportionally.  However, 
when  subsample  size  differences  are  relatively  small  so  that 
overlap  of  the  box  plots  would  occur,  a  minimum  of  six  spaces 
between  centers  of  box  plots  is  left  to  prevent  the  graphics 
from  becoming  confused. 

The  vertical  scale  of  the  graph  is  determined  by  the 
minimum  and  maximum  values  from  the  population  of  estimates 
determined  from  the  smallest  subsample  size  (n^).  It  is 
assumed  here  that  subsequent  subsample  sizes  (n2 ,  n-,  ...) 
will  be  large  in  relation  to  the  smallest  subsample  size 
Cn^^)  so  that  estimates  produced  from  these  will  be  within 
the  range  of  that  determined  from  the  smallest  subsample 
size.  This  assumption  appears  tenable  if  the  range  of  esti¬ 
mates  is  compared  in  Figure  4a  for  increasing  subsample  sizes. 
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The  vertical  scale  of  Figure  4a  is  not  scaled  to  the 
minimum  and  maximum  value  of  the  estimates  determined  from 
the  smallest  subsample  size  (nj^-33) .  This  is  because  the 
RAGE  program  can  evaluate  up  to  three  different  estimators 
at  one  time  and  two  alternate  estimators  (Figures  4b  S  4c) 
were  also  used  to  estimate  the  parameter.  An  option  was 
entered  to  scale  all  three  estimators  to  the  same  scale  and 
it  chose  the  minimum  and  maximum  values  from  all  the  esti¬ 
mates  produced  by  all  three  estimators  using  the  smallest 
subsample  size  (n2^*33).  By  comparing  Figures  4a,  4b  and  4c 
it  can  be  seen  that  the  maximum  value  was  taken  from  the 
first  estimator,  Figure  4a,  and  the  minimum  value  was  taken 
from  the  third  estimator.  Figure  4c.  Other  scaling  options 
will  be  discussed  in  detail  in  the  following  chapter. 

Below  the  graph  statistical  results  as  they  pertain  to 
the  population  of  estimates  for  each  subsample  size  are 
printed.  These  statistics  include  the  sample  mean,  standard 
deviation,  standard  deviation  of  the  mean,  skewness  and  kur- 
tosis.  Thus,  for  example,  in  Figure  4a  at  subsample  size 
n^“33  there  is  a  population  of  [N/nj^]xM  »  rj^xM  »  [500/33]xl00  » 
150  estimates  for  which  a  sample  mean,  standard  deviation, 
standard  deviation  of  the  mean,  skewness  and  kurtosis  is 
computed.  The  subsample  sizes  these  statistics  refer  to  are 
labeled  across  the  top. 

As  mentioned  earlier,  in  addition  to  box  plots  and  associ¬ 
ated  statistics  of  the  estimates,  RARE  and  unbiased  RARE 
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estimates  of  the  parameter  are  computed  and  plotted  on  the 
graph.  The  degree  of  regression,  D,  for  the  computations  of 
the  unbiased  RARE  estimates  is  determined  by  the  user  and 
printed  at  the  top  of  the  graph  for  future  reference.  The 
RARE  estimate  is  the  regression  coefficient  3q  and  is  the 
estimated  asymptotic  expected  value  of  the  parameter  being 
estimated.  The  unbiased  RARE  estimates  are  determined  by 
the  regression  coefficients,  $2,  ...,  the  subsample 

size,  n,  and  the  RARE  estimate,  Sg,  using  the  following 
formula: 

Sg  >  32-^  ^  .  • 

n  “  n° 

A  set  of  regression  coefficients  (3g,  8^,  ....  3^)  is 
determined  from  the  mean  of  the  population  or  estimates 
(dependent  variable)  at  each  subsample  size  for  each  replica 
tion.  The  M  regressions  that  determine  the  M  sets  of  coeffi 
cients  (Sg.  ...,  S^j)  are  represented  by  the  following 

formula: 

m^(ni)  -  8o  *  ^  .  ^  A 

1  nr  n“ 

1  1 


m  ■  1 ,  . . . ,  M 
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K 

M 


A 

where  S  (n,  ) 
m  k 


^  A 


k  =  1 ,  . . . , 

m  “  I ,  •  •  •  > 


The  regression  coefficients,  to  determine  the  RARE  and  unbiased 
RARE  estimates,  are  the  average  values  determined  from  these 
M  regressions  applied  to  the  M  replications.  The  mean, 
variance  and  standard  deviation  of  these  coefficients  are 
printed  below  the  graph  below  the  subsample  size  statistics. 

As  mentioned  earlier,  the  first  coefficient,  3o,  is  the 
estimated  asymptotic  expected  value  of  the  parameter  being 
estimated  (RARE).  Its  mean  from  the  M  regressions  is  printed 
as  a  horizontal  dashed  line  across  the  graph.  As  the  sub- 
sample  size  increases  the  bias  in  the  unbiased  RARE  estimates 
becomes  smaller  and  the  unbiased  RARE  estimates  (dotted  line) 
will  asymptotically  approach  the  RARE  estimate  (dashed  line) . 

In  Figure  4a  the  estimated  expected  value  or  asymptote 

Sq  *  4.98  and  the  bias  coefficients  for  the  three  degree  (D) 

regression  are:  =»  16.50,  83  “  -182.21,  63  *  3659.04. 

The  variance  regression  coefficients  are  determined  from 
the  variances  of  the  population  of  estimates  for  each  sub¬ 
sample  size  (dependent  variable)  and  evaluating  the  single 
regression 
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The  four  (D+1)  variance  regression  coefficients  for  Figure  4a 
are;  aj^  •  58.37,  a2  =  174.72,  =  -1500.13  ,  and  a^  *  5446.25. 

The  variance  regression  coefficients  derived  here  can  be  used 
to  obtain  an  approximate  variance  for  a  specified  subsample 
size  by  substituting  the  subsample  size  for  nj^  in  any  one 
of  the  above  regression  equations  along  with  the  variance 
coefficients  (a^,  02,  ..•)  and  completing  the  computation. 

Below  the  variance  regression  coefficients  the  minimum 
and  maximum  value  of  the  vertical  scale  is  printed.  The 
purpose  of  this  information  is  to  allow  the  user  to  scale 
other  estimators  of  the  same  parameter  to  the  same  scale. 

It  can  be  used  when  the  number  of  estimators  exceeds  three 
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(the  maximum  number  of  estimators  RAGE  can  use  in  a  single 
run)  or  when  subsequent  runs  for  different  conditions  are 
desired  to  be  of  the  same  scale  for  graphical  comparisons. 

A  final  identifying  feature  of  the  graphics  is  a  section 
at  the  very  bottom  of  the  graphical  output  which  allows  the 
user  to  identify  the  estimator  that  was  used  in  estimating 
the  parameter.  In  Figure  4a  this  estimator  was  the  moment 
estimator,  i.e.  the  reciprocal  of  the  estimated  coefficient 
of  variation  squared. 
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III.  DESCRIPTION  OF  RAGE  ARGUMENTS 


The  RAGE  program  has  the  calling  sequence:  call 
RAGE(X,N,M,NE,L,D,RG,SEI ,SVS ,YMIN,YMAX,NEST,EST1 ,TTL1 ,EST2 , 
TTL2 ,EST3 ,TTL3) .  The  arguments  are  relatively  easy  to 
understand  but  in  several  cases  are  related  and  must  meet 
certain  restrictions.  For  this  reason  a  detailed  descrip¬ 
tion  of  each  argument  is  given,  as  well  as  any  restrictions 
or  limitations  that  must  be  met  within  or  between  arguments. 
Following  this  detailed  discussion  an  example  is  given  in 
Chapter  IV  on  how  to  set  up  the  main  program  and  the  estima¬ 
tor  functions. 

A.  DATA  ARRAY. . .X 

The  X  argument  is  a  single  precision  Real*4  array  con¬ 
taining  the  data  from  which  the  parameter  is  to  be  estimated. 
The  only  restriction  is  that  the  size  of  the  X  array  must 
not  exceed  50,000  which  is  the  maximum  sample  size  that  RAGE 
can  store. 

B.  SAMPLE  SIZE. . .N 

The  number  of  data  elements  per  section  or  sample  size 
is  designated  by  the  integer  N.  Depending  on  the  minimum 
sample  size  required  by  the  estimator,  N  can  vary  from  1  to 
50,000.  The  sample  size,  N,  times  the  number  of  sections, 

M,  cannot  exceed  50,000  (the  maximum  number  of  data  elements 
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stored  in  the  array  X).  If  M  times  N  exceeds  50,000  an 
error  message  will  be  printed  and  execution  will  terminate. 

If  M  times  N  exceeds  the  number  of  data  elements  stored  by 
the  user  in  the  array  X  then  no  error  message  will  be  printed 
and  the  program  will  execute  giving  unpredictable  results  or 
an  abnormal  ending  will  occur. 

C.  NUMBER  OF  REPLICATIONS .. .M 

The  number  of  sections  or  replications  the  array  X  is 
sectioned  into  is  determined  by  the  integer  value  assigned  to 
M.  The  number  assigned  to  M  is  also  the  number  of  regressions 
that  will  be  run  to  determine  the  regression  coefficients 
Sq,  ..Sji.  The  number  of  regressions  is  limited  to  one 

hundred  and  therefore  M  is  limited  to  a  maximum  value  of  one 
hundred.  An  error  message  will  be  printed  and  execution  will 
terminate  if  M  is  greater  than  one  hundred  or  less  than  one. 

If  M  is  equal  to  one  the  variance  and  standard  deviation  of 
the  regression  coefficients  for  the  unbiased  estimate  will 
not  be  printed  since  they  will  not  exist. 

D.  SUBSAMPLE  SIZE  ARRAY... NE 

The  argument  NE  is  an  integer  array  of  size  eight  con¬ 
taining  the  subsample  sizes  (n^^,  k=l . K)  the  estimator 

is  to  use  in  estimating  the  parameter.  The  values  in  the 
array  NE  should  be  ordered  from  smallest  to  largest  or  a 
warning  message  will  be  issued  to  indicate  the  vertical 
scale  of  the  graph  may  not  be  large  enough  to  include  all 
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estimates.  RAGE  will  use  the  population  of  estimates  deter¬ 
mined  from  the  subsample  size  stored  in  NE(1)  to  determine 
the  vertical  scale  of  the  graph.  If  NE(1)  is  not  the  smallest 
value  the  range  of  the  graphics  may  not  be  large  enough  to 
include  outliers  determined  from  subsequent  smaller  subsample 
sizes . 

Due  to  storage  limitations  a  maximum  of  12,500  estimates 
is  perraissable.  Since  NEfl)  will  produce  the  most  estimates, 
the  user  can  determine  the  total  number  of  estimates  using  the 
formula  Mx[N/NE(l)].  If  this  value  exceeds  12,500  an  error 
message  will  result  stating  this  limitation  has  been  exceeded 
and  execution  will  terminate.  Corrective  action  can  be  taken 
by  increasing  NE(1)  or  decreasing  M  or  N. 

E.  NUMBER  OF  SUBSAMPLE  SIZES  (BOX  PLOTS)... L 

The  user  can  select  the  number  of  subsample  sizes  to 
be  used  from  the  array  NEC8)  by  assigning  an  integer  value 
to  the  argument  L.  The  argument,  L,  can  take  on  values 
between  one  and  eight.  RAGE  will  then  select  the  first  L 
subsample  sizes  from  NE(8),  compute  the  population  of  esti¬ 
mates  of  the  parameter  with  each,  graph  the  box  plots  and 
print  the  associated  statistics.  If  L  is  not  an  integer 
value  between  one  and  eight  an  error  message  will  result  and 
execution  will  terminate. 

If  L  =  1  no  regressions  are  possible  and  any  regression 
output  is  omitted.  The  maximum  degree  of  regression  allowed 
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by  RAGE,  regardless  of  the  value  assigned  to  the  argument  D 
(degree  of  regression),  is  L-1.  If  D  exceeds  this  value  RAGE 
will  set  the  degree  of  regression  to  L-1. 

F.  DEGREE  OF  REGRESSION. . .D 

The  degree  of  regression,  D,  is  an  integer  value  less 
than  or  equal  to  three.  The  degree  of  regression  refers  to 
the  number  of  bias  coefficients,  3,,  ..3^,  that  will  be 

used  in  computing  the  unbiased  estimate  out  to  D  terms. 

In  addition  to  the  restrictions  on  D  mentioned  earlier, 

D  will  also  be  reduced  for  the  regression  that  determines 
the  variance  coefficients,  a^,  cx^,  when  the  number  of 

subsamples  sites  with  a  variance  is  less  than  D+1. 

G.  REDUCED  GRAPHICS. . .RG 

The  reduced  graphics  argument,  RG,  allows  the  user  to 
reduce  the  range  of  the  vertical  scale  of  the  graph  by  elimi¬ 
nating  outliers  and  only  printing  the  number  of  outliers 
that  exist.  The  vertical  scale  is  reduced  to  the  upper  quar- 
tile  plus  1.5  times  the  interquartile  distance  value  and  the 
lower  quartile  minus  1,5  times  the  interquartile  distance 
value  of  the  estimates  from  the  smallest  subsample  size, 
NE(1).  The  vertical  scale  then  becomes  these  values  and  the 
outliers  are  merely  counted  and  printed  on  the  graph  at  the 
ends  of  each  box  plot's  upper  (lower)  quartile  +  (-)  1.5 
times  interquartile  distance  position.  If  there  are  no  out¬ 
liers  for  the  NE(1)  estimates,  either  or  both  extremes,  the 
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vertical  scale  is  further  reduced  to  the  first  estimate 
value  within  these  limits.  RG  »  1  will  cause  reduced  graphics 
to  be  produced  and  RG  =  0  will  cause  full  scale  graphics. 

H.  SCALING  ESTIMATORS  INDIVIDUALLY ... SEI 

The  argument  SEI  allows  the  user  to  scale  each  estimator's 
vertical  scale  for  the  graphics  separately  or  to  the  same 
scale.  SEI  =  0  causes  all  estimators'  vertical  scale  to  have 
the  same  range  of  values.  The  range  of  the  vertical  scale 
and  whether  outliers  are  printed  or  counted  will  depend  on  the 
graphics  desired  (RG  =  0  or  1).  Figures  4a,  4b  and  4c  are  the 
results  of  three  estimators  scaled  to  the  same  scale  (SEI  =  0) 
with  full  scale  graphics  (RG  *  0).  Figures  5a,  5b  and  5c  are 
the  results  of  three  estimators  scaled  to  the  same  scale 
(SEI  »  0)  but  with  reduced  graphics  (RG  *  1). 

If  SEI  =»  1  each  estimator's  vertical  scale  will  be  scaled 
individually.  The  range  of  the  vertical  scale  and  whether 
outliers  are  printed  or  counted  will  depend  on  the  graphics 
desired  (RG  *  0  or  1) . 

I.  SETTING  THE  VERTICAL  SCALE... SVS 

To  set  the  vertical  scale  an  integer  value  of  1  is  assigned 
to  the  argument  SVS.  SVS  =  1  will  cause  the  SEI  argument  to 
be  ignored  and  the  vertical  scale  of  all  the  estimators' 
graphics  will  be  set  to  those  values  assigned  to  YMIN  and 
YMAX.  When  the  argument  SVS  is  in  effect  (SVS  *  1)  the 
argument  RG  is  also  in  effect  as  far  as  whether  outliers 
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are  printed  or  counted  (RG  »  0  or  1)  but  the  determination 
of  the  vertical  scale  by  RG  is  overridden  by  SVS  and  is  set 
to  YMIN  and  YMAX.  SVS  *  0  causes  RAGE  to  ignore  values 
assigned  to  YMIN  and  YMAX  and  scaling  is  determined  by  the 
values  assigned  to  the  arguments  RG  and  SEI.  Figure  4d  is 
an  example  of  SVS  *  1  and  RG  »  0  with  the  vertical  scale 
set  to  that  of  Figures  4a,  4b,  and  4c  (YMIN=-1 . 2812 , 

YMAX*13 . 9984) .  Figure  Sd  is  an  example  of  SVS  »  1  and 
RG  »  1  with  the  vertical  scale  set  to  that  of  Figures  5a, 

5b  and  5c  (YMIN=1 . 0394 ,  YMAX=9. 1423) . 

J.  MINIMUM  VALUE  OF  THE  VERTICAL  SCALE... YMIN 

The  value  YMIN  is  a  real  number  assigned  by  the  user  as 
the  lower  limit  of  the  vertical  scale.  If  YMIN  is  assigned 
a  value  that  is  so  large  that  it  falls  within  the  range  of 
values  that  determine  the  body  of  any  box  plot  then  an  error 
will  occur  and  an  abnormal  ending  will  result. 

K.  MAXIMUM  VALUE  OF  THE  VERTICAL  SCALE... YMAX 

The  value  YMAX  is  a  real  number  assigned  by  the  user  as 
the  upper  limit  of  the  vertical  scale.  If  YMAX  is  assigned 
a  value  that  is  so  small  that  it  falls  within  the  range  of 
values  that  determine  the  body  of  any  box  plot  then  an  error 
will  occur  and  an  abnormal  ending  will  result. 

L.  NUMBER  OF  ESTIMATORS .. .NEST 

The  argument  NEST  is  an  integer,  1,  2  or  3,  which  indicates 
the  number  of  estimators  that  will  be  used  to  estimate  the 
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parameter.  An  error  message  will  be  issued  if  NEST  is  not 
one  of  these  values  and  execution  will  terminate. 

M.  ESTIMATOR  FUNCTIONS. . .ESTl ,EST2 ,EST3 

These  arguments  are  the  function  names  assigned  to  each 
of  the  estimators  that  are  to  be  used  in  estimating  the  para 
meter.  Each  function  must  be  written  with  a  two  argument 
calling  sequence  (i.e.  call  EST1(X,N))  where  the  first 
argument,  X,  must  be  the  data  array  and  the  second  argument, 

N,  must  be  the  number  of  data  elements. 

If  there  are  less  than  three  estimators  then  dummy  argu¬ 
ments  must  be  inserted  in  the  calling  sequence  where  subse¬ 
quent  estimators  would  normally  appear  (Suggestion. . .use  a 
function  of  a  previous  estimator*,  i.e.  ESTl,  as  the  dummy. 
RAGE  will  not  provide  duplicate  output  for  these  repeated 
functions  unless  NEST  is  increased  to  include  them.).  Each 
function  must  then  be  declared  external  in  the  main  program. 
For  an  illustration  of  how  the  functions  are  written  and 
declared  see  the  example  in  Chapter  IV. 

N.  DESCRIPTION  OF  ESTIMATORS. . .TTLl ,TTL2 ,TTL3 

The  titles  associated  with  each  estimator  can  be  passed 
as  an  argument  through  the  call  program  (RAGE)  provided 
there  is  a  minimum  of  120  characters  between  apostrophes 

(i.e.  call  RAGE(...,  ESTl,' - 120  characters - ',EST2,...)) 

An  alternate  and  possibly  easier  method  is  to  declare  each 
title  (TTLl,  TTL2,  TTL3)  as  a  Real*8(15)  array  and  assign 
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the  characters  that  describe  each  estimator  as  data  to  each 
title.  For  an  example  of  this  method  see  the  example  in 
Chapter  IV. 

The  user  should  note  that  the  title  associated  with  each 
estimator  must  immediately  follow  the  estimator  argument. 

As  in  the  case  of  estimators,  dummy  titles  must  be  inserted 
when  the  number  of  estimators  requiring  titles  is  less  than 
three.  Finally,  it  should  be  noted  that  the  program  will 
not  execute  properly  unless  titles  of  the  correct  format  are 
included. 
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IV.  EXAMPLE  OF  A  PROGRAM  USING  RAGE 


In  the  following  example  suppose  a  user  has  a  set  of  iid 
data  from  which  he  wishes  to  estimate  the  population  mean  by 
using  two  estimators,  the  sample  mean  and  the  sample  median. 
Assume  that  these  estimators  are  given  the  abbreviated  func¬ 
tion  names  MN  and  MED  and  that  the  associated  titles  for 
each  estimator  are  "The  sample  mean  as  an  estimator  of  the 
population  mean"  and  "The  sample  median  as  an  estimator  of 
the  population  mean."  If  the  subsample  sizes  at  which  it  is 
desired  to  examine  the  properties  of  these  estimators  are 
33,  so,  71,  12S,  166,  200,  250  and  500  the  FORTRAN  program 
may  look  as  formatted  on  the  following  page. 
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REAL*4  XCSOOOO) .YMIN.YMAX 
REAL*8  T1 (15) .12(15) 

INTEGER  NEST,D,RG,SEI,SVS,M.N,L,NE(8) 

DATA  NE/  33,50,71,125,166,200,250  500/ 

C  INSERT  8  CHARACTERS  BETWEEN  EACH  SET  OF  APOSTROPHES  UNTIL 

C  TITLE  IS  COMPLETED  AND  FILL  REMAINDER  WITH  BLANKS. 

DATA  Tl/'THE  SAMP ' , ’ LE  MEAN  ','AS  AN  ES'.'TIMATOR  ', 

+'OF  THE  P' ,'OPULATIO' ,'N  MEAN', 8*'  '/ 

DATA  T2/'THE  SAMP'.'LE  MEDIA’, 'N  AS  AN  ' , ' ESTIMATO ' , 

♦’R  OF  THE’ ,'POPULATI' ,’ON  MEAN’, 8*’  ’/ 

EXTERNAL  MN.MED 


CALL  RAGE(X,N,M,NE,L,D,RG,SEI ,SVS ,YMIN , YMAX , 
+NEST,MN,T1 ,M5D,T2 ,MN,T1) 


STOP 

END 
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REAL  FUNCTION  MN(X,N) 

REAL  X(N) 

C  AVERAGE  OF  N  VALUES  IN  X 

MN-0.0 
DO  10  I-1,N 

MN  •  MN  +XCI) 

10  CONTINUE 

MN  -  MN/ FLOAT (N) 

RETURN 

END 

REAL  FUNCTION  MED(X,N) 

REAL  X(N) 

C  MEDIAN  OF  N  VALUES  IN  X 

C  "PCTILE”  IS  A  FUNCTION  THAT  RETURNS  THE  SO  PERCENTILE 

C  NUMBER  FROM  THE  DATA  ARRAY  X  OF  SIZE  N. 

MED  •  PCTILE(X,N,.5) 

RETURN 

END 

In  this  example  since  only  two  estimators  were  used  a  dummy 
estimator  and  title  had  to  be  inserted  as  the  arguments  for 
the  third  case.  Here  the  first  estimator  and  its  title  were 
simply  repeated. 
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r.  A  COMPARATIVE  STUDY  OF  ESTIMATORS  FOR  THE  GAMMA  SHAPE 

PARAMETER 

To  demonstrate  the  use  o£  the  Regression  Adjusted  Graphics 
and  Estimating  (RAGE)  package,  two  arrays  of  gamma  distri¬ 
buted  random  variables  with  shape  K»5.0  and  K»0.2  were 
generated  using  the  Naval  Postgraduate  School's  random  number 
generator- -LLRANDOMI I  [Ref.  3].  Because  these  random  numbers 
are  independent  and  were  drawn  from  the  same  distribution, 
they  can  be  sectioned  into  independent  and  identically 
distributed  blocks  of  size  N.  RAGE  was  specifically  designed 
to  evaluate  estimators  of  parameters  from  just  such  iid 
populations.  These  two  samples  then  should  demonstrate  the 
capabilities  of  the  RAGE  package  to  evaluate  the  bias,  vari¬ 
ance  and  distribution  of  estimators  of  a  specified  parameter 
from  the  populations. 

In  the  following  examples  it  is  the  shape  parameter,  K, 
of  the  gamma  distribution  that  is  being  estimated.  There  are 
two  primary  estimators  used  in  these  examples.  The  first  is 
the  moment  estimator  (reciprocal  of  the  estimated  coefficient 

1  Y- 

of  variation  squared),  K  =  _  =  _ _  .  The  second 

(C(X))2  VAR(X) 

estimator  is  the  maximum  likelihood  estimate  (MLE)  [Refs.  2,4] 
which  is  determined  by  a  bisection  search  for  the  shape  para¬ 
meter,  K,  until  the  following  equation  is  satisfied: 
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Ln(k)  -  Diganuna(K)  »  Ln  X  - 


1-1 


n 


i=l ,  .  .  . , 


n . 


The  i*l,  2 . n  are  the  random  numbers  from  a  specified 

block  of  size  n.  The  other  two  estimators  investigated  are 
jacknifed  versions  of  the  initial  estimators.  The  jacknife 
routine  [Ref.  5]  with  each  of  these  estimators  is  a  four  fold 
jacknife  which  has  the  form 


’w 

0 


* 


1 

4 


S* 


j 


"Vi  'X,  a,  'X, 

where  -  ^(®all^  '  ®all  estimate  from  the 

-v. 

entire  sample  (n)  and  0j  is  the  estimate  with  the  j^j^  group 
missing.  Each  group  is  a  set  of  one  fourth  of  the  data 
points . 

Each  of  the  four  estimators  is  examined  for  varying  sub¬ 
sample  sizes  and  replications  to  determine  how  accurately 
each  estimator  estimates  the  shape  parameter  of  the  gamma 
distribution  as  the  subsaraple  size  increases.  The  shape 
parameter  of  the  gamma  distribution  was  chosen  as  an  example 
because  the  behavior  of  the  many  estimators  which  have  been 
proposed  for  K  is  fairly  well  known.  Roughly  speaking,  it  is 
known  that  the  moment  estimator  is  quite  efficient  for  K>2 
but  quite  inefficient  for  K<1  even  for  large  sample  sizes. 

Bias  and  normality  of  these  estimators  is  not  documented. 

Again  the  jacknife  is  known  to  work  well  with  moment  estimators 


even  in  small  samples  but  no  such  results  are  known  for  the 
jacknifed  MLE  estimator  except  that  it  is  asymptotically 
efficient . 


For  this  investigation  there  are  basically  three  examples 
representing  different  simulation  scenarios  for  evaluating 
each  of  two  cases,  K»5.0  and  K*0.2.  They  are: 

(a)  Small  subsample  sizes  (33  to  500)  ,  small 
replications  (M=2) ; 

(b)  Small  subsample  sizes  (33  to  500)  ,  large 
replications  (M=100) ; 

(c)  Large  subsample  sizes  (100  to  2000)  ,  large 
replications  (M=2S). 

The  first  can  be  looked  on  as  a  pilot  study  to  determine 
regions  of  interest  for  the  two  longer  scale  studies.  How¬ 
ever,  the  procedures  are  not  fully  sequential,  unlike  the 
usage  of  the  RAGE  methodology  proposed  by  Heidelberger  and 
Lewis  [Ref.  1]  for  regenerative  simulations. 

With  many  of  these  examples  another  set  of  reduced  graphics 
is  included  which  reduces  the  vertical  scale  and  provides 
graphics  which  are  more  suitable  for  examining  the  bias 
structure  of  the  estimators.  The  reduced  graphics  truncates 
the  outliers  in  the  box  plots  giving  only  their  number  and 
not  their  position. 
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A.  EXAMINING  THE  GAMMA  DISTRIBUTION  FOR  THE  SHAPE  PARAMETER 
(K=5.0) 

1 .  Case  I.  Example  1:  Small  Sub sample  Sizes  and 

Small  Replications  CK=5.01 

Figures  3a,  3b,  3c  and  3d  show  the  graphical  and 
statistical  results  from  estimating  the  shape  parameter  with 
the  four  different  estimators.  These  graphs  all  have  the 
same  vertical  scale  so  that  meaningful  comparisons  can  be 
made  of  the  box  plots  and  of  the  bias  of  each  estimator. 

One  of  the  characteristics  of  the  MLE  estimator  is 


that  its  estimates  will  asymptotically  have  the  smallest  vari¬ 
ance  of  all  estimators.  For  small  subsample  sizes  (such  as 
the  case  of  Figures  3a  §  3b)  the  asymptotically  smaller  vari¬ 
ance  of  the  MLE  estimator  is  not  apparent  in  the  box  plots. 

The  body  of  the  box  plot  for  the  MLE  estimator  for  n^  =  33 
is  smaller  than  the  same  box  plot  for  the  moment  estimator. 
However,  the  MLE  has  more  and  larger  outliers,  indicating 
that  the  MLE  is  more  nonnormal  at  this  sample  size  than  the 
moment  estimator.  If,  however,  the  estimated  standard  devia¬ 
tion  of  the  various  subsample  sizes  are  compared  for  each  of 
the  two  primary  estimators  (Table  1)  it  can  be  seen  that  the 
MLE  estimator  consistently  produces  a  smaller  variance  than 
does  the  moment  estimator.  From  Figures  3a  5  3b  the  bias  in 
the  two  estimators  appears  comparable,  starting  high  but 
becoming  negligible  by  subsample  size  250. 


Table  1;  Standard  deviations  of  estimates  of  the  shape 
parameter  using  the  moment  and  MLE  estimators 
for  small  sample  sizes  (replications  =  2) 


Standard 

Deviation 

SAMPLE  SIZE 

MOME.MT  ESTIMATOR 

MLE  ESTIMATOR 

33 

1.517 

1 . 306 

so 

.9942 

.7151 

71 

.  8829 

.  6051 

100 

.6312 

.4612 

125 

.  7725 

.4823 

166 

.8310 

.5114 

250 

.4660 

.  2772 

500 

.5589 

.  3087 

Jacknifing  is  a  technique  that  can  be  used  on  iid  sam¬ 
ples  to  reduce  the  bias  in  the  estimates  caused  by  the  estimator 
[Ref.  5].  An  estimator  with  a  complete  jacknife  is  asymptoti¬ 
cally  normal  with  the  same  variance  as  the  unjacknifed  estimator, 
but  small  sample  properties  are  difficult  to  obtain  analyti¬ 
cally.  Figures  3c  §  3d  are  the  jacknifed  estimators  and 
again  because  of  the  small  number  of  estimates  it  is  not 
apparent  that  the  jacknife  is  contributing  to  nonnormality 
or  increasing  variance  when  the  box  plots  are  compared  with 
those  of  Figures  3a  §  3b.  The  increased  normality,  however, 
can  be  seen  in  the  jacknifed  moment  estimator  by  comparing 
box  plots  (Figures  3a  §  3c)  and  the  coefficients  of  skewness 
(Table  2)  with  the  nonjacknifed  form.  The  change  in  normality 
or  symmetry  of  the  estimates  from  the  jacknifed  MLE  estimator 
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is  not  discernable  in  a  comparison  of  the  box  plots  (Figures 
3b  §  3d)  or  the  coefficients  of  skewness  (Table  3). 

Table  2:  Skewness  coefficients  of  the  estimates  of 
the  shape  parameter  using  the  moment  and 
jacknifed  moment  estimators  (replications  =  2) 

Skewness  Coefficients 

SAMPLE  SIZE  MOMENT  ESTIMATOR  JACKNIFED  FORM 


33 

.553 

.327 

50 

.076 

-.348 

71 

.339 

.253 

100 

.499 

.157 

125 

1.573 

.724 

166 

1.018 

.501 

250 

-.150 

-  .235 

500 

N/A 

N/A 

Table  3:  Skeimess  coefficients  of  the  estimates  of  the 

shape  parameter  using  the  MLE  and  jacknifed  MLE 
estimators  (replications  =  2) 

Skewness  Coefficients 

SAMPLE  SIZE  MLE  ESTIMATOR  JACKNIFED  FORM 


33 

1.198 

.050 

50 

.510 

.616 

71 

.193 

-.163 

100 

.773 

.  643 

125 

1.451 

1.216 

166 

1.101 

1.229 

250 

-.030 

.081 

500 

N/A 

N/A 
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The  bias  in  each  estimator  is  depicted  by  the  space 
between  the  dotted  line  and  dashed  line  of  each  of  the  graphs. 
When  the  unbiased  estimate  Cd-otted  line)  approaches  the 
estimated  e.\pected  value  (dashed  line)  then  the  bias  becomes 
insignificant.  Each  estimator  can  be  expected  to  vary  in 
the  amount  of  bias  it  has  in  its  estimates  and  this  bias  can 
be  expected  to  diminish  as  sample  size  increases.  Comparing 
Figures  3a  §  3b  it  appears  that  the  MLE  estimator  has  about 
the  same  bias  as  the  moment  estimator  but  this  is  hard  to 
determine  at  this  point  since  the  bias  coefficients  were 
determined  from  a  small  number  of  estimates.  In  the  esti¬ 
mates  from  the  jacknifed  estimators  one  would  expect  the  bias 
to  be  reduced  because  jacknifing  is  specifically  designed  to 
remove  the  (1/n)  term  in  the  bias  of  the  estimates  but 
this  is  not  obvious  from  the  graphical  comparisons  (Figures 
3a  S  3c  and  Figures  3b  S  3d).  For  the  case  of  the  MLE  esti¬ 
mator  the  coefficient  of  the  1/n  term,  ,  has  actually 
increased  from  32.65  to  75.33  after  jacknifing.  The  standard 
deviation  of  these  bias  terms  are  so  large  (19.14  and  16.23 
respectively),  though,  that  they  suggest  more  estimates 
(increased  replications)  are  necessary  to  better  determine 
the  bias. 

Determining  how  much  more  bias  one  estimator  has  than 
another  can  also  be  accomplished  by  comparing  the  bias 
regression  coefficients  Bj^,  62#  etc.  The  larger  these  coef¬ 
ficients  are  the  greater  is  the  bias  of  that  particular 
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estimator.  Comparing  the  coefficients  for  the  moment  and 
MLE  estimators  of  Figures  3a  §  3b  show  that  the  MLE  estimator 
has  a  larger  bias.  For  a  subsample  size  of  500  and  a  three 


degree 

regression 

the  bias  for 

the  moment 

estimator 

7.40128 

183.099 
+  + 

1947.86 

=  .0155 

500 

500^ 

500^ 

and  the 

bias  for 

the  MLE  estimator  is 

32.6486 

1998.68 

47308.0 

=  .0577 

500 

500^ 

500^ 

These  comparisons  seem  to  indicate  that  the  MLE  estimator 
has  more  bias.  This  is  not  necessarily  true,  though,  for 
reasons  mentioned  earlier.  Increasing  the  number  of  replica¬ 
tions  will  alleviate  much  of  this  uncertainty  since  more 
precision  in  the  estimates  of  the  bias  coefficients  will 
result . 
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SAHPli  SUt  INI:  »00  NU.  Uf  KEPLlCAlIUNS  INI:  2  1>E6R£E  Of  REMEE&IUN  10|t 


2.  Case  1,  Example  2:  Small  Subsample  Sizes  and  Large 


It  was  pointed  out  in  Example  1  of  this  case  that 
the  small  number  of  replications  left  some  doubt  as  to  the 
accuracy  of  the  estimates  of  the  bias  coefficients  as  well 
as  the  normality  of  the  estimates.  Figures  4a  to  4d  (cor¬ 
responding  to  Figures  3a  through  3d)  are  for  increased 
replications  (M=100).  This  has  the  effect  of  increasing  the 
population  of  estimates  for  each  subsample  size  (nj.,  k»l,  .... 
K)  and  the  number  of  estimates  of  the  bias  coefficients 
(3i»  62*  holds  the  subsample  size  at  which  the  esti¬ 

mators  are  examined  small  as  before.  On  first  observing  the 
graphs,  it  becomes  apparent  from  the  box  plots  of  the  two 
smallest  subsample  sizes  Cnj^*33  and  n2*50)  that  the  estimates 
are  nonnormally  distributed  (Figures  4a  5  4b).  The  jacknifed 
estimators  remove  some  of  this  nonnormality  (Figures  4c  §  4c) , 
especially  for  the  jacknifed  moment  estimator.  In  effect, 
outliers  appear  now  at  both  high  and  low  values  but  consider¬ 
able  skewness  still  remains. 

Here  it  is  apparent  by  the  statistical  output  that  the 
MLE  estimator  has  significantly  less  variance  in  its  popula¬ 
tion  of  estimates  than  does  the  moment  estimator  (Table  4). 


Table  4;  Standard  deviations  of  the  estimates  of  the 
shape  parameter  using  the  moment  and  MLE 
estimators  for  varying  sample  sizes 
(replications  »  100) 


Standard 

Deviation 

SAMPLE  SIZE 

MOMENT  ESTIMATOR 

MLE  ESTIMATOR 

33 

1.486 

1.396 

SO 

1.175 

1.075 

71 

.9818 

.  8549 

100 

.8008 

.6636 

12S 

.7242 

.5918 

166 

.6499 

.5145 

2S0 

.  5082 

.4036 

500 

.3417 

.2730 

The  increased  symmetry  and  normality  by  the  jacknifed  moment 
estimator  is  graphically  visible  if  Figures  4a  §  4c  are  com¬ 
pared.  The  change  in  symmetry  or  normality  of  the  estimates 
of  the  jacknifed  MLE  estimator  was  not  discernable  in  the 
graphics  (Figures  3b  5  3d)  or  skewness  coefficients  (Table  3) 
for  Example  1.  However,  for  increased  replications,  (M*100) 
if  the  box  plots  of  Figures  4b  §  4d  are  compared  it  appears 
that  the  population  of  MLE  estimates  are  more  normally 
distributed  (outliers  at  both  ends)  after  jacknifing.  The 
movement  toward  symmetry  is  more  obvious  by  comparing  the 
skewness  coefficients  of  the  MLE  estimates  and  the  jacknifed 
MLE  estimates. 
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The  variance  reduction,  or  lack  of  it,  for  each 
estimator  when  it  is  jacknifed  can  be  determined  by  comparing 
the  computed  variances  that  result  from  a  given  subsample 
size  and  the  variance  coefficients,  a2 ,  etc.  If,  for 
example,  the  variance  coefficients  (^1.  a2 »  '^3» 
each  of  the  estimators  (Figures  4a  through  4d)  are  used  to 
determine  the  variance  for  a  subsample  size  of  500,  the  vari¬ 
ances  for  each  estimator  would  be: 


VAR (MOMENT)  » 

58.37 

+ 

174.72 

1500.13 

5446.25 

+  a 

.1273 

500 

500^*^ 

500^ 

500^ • ^ 

VAR (MLE) 

67.52 

961.91 

10183.5 

27503.3  = 

.  0848 

500 

500^*^ 

500^ 

500^*^ 

VAR(JNK  MNT)- 

79.20 

457.44 

+ 

5216.8 

14533.1 

.1357 

500 

500^*5 

500- 

500^  *  ^ 

VAR(JNK  MLE)» 

62.66 

825.2 

+ 

8879.79 

24484.2 

3 

.  0827 

500 

500^-^ 

500^ 

500^*^ 

These  computations  show 

a  lack 

of 

change 

in  the  variance 

when  the  estimators  are  jacknifed  and  agree  with  the  statis¬ 
tical  and  graphical  conclusions  drawn  earlier. 

The  bias  in  the  MLE  and  moment  estimators  is  difficult 
to  compare  graphically  (Figures  4a  §  4b)  since  the  unbiased 
estimate  approaches  the  estimated  expected  value  so  quickly 
and  because  of  the  compacting  caused  by  the  range  of  the 
vertical  scale.  For  this  reason  reduced  graphics  are  produced 
(Figures  5a  through  5d) ,  which  the  reader  will  recall 
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suppresses  the  outliers.  Since  the  scale  is  reduced  from  a 
range  of  15,28  (14.00  -  (-1.28))  to  8.10  (9.14  -  1.04)  or 
almost  50%,  the  evolution  of  the  bias  as  subsample  size 
increases  becomes  clearer,  as  Figures  Sa  8  5b  demonstrate  in 
comparison  to  Figures  4a  S  4b.  These  graphs  indicate  that, 
for  at  least  small  subsample  sizes,  the  moment  estimator  is 
a  more  biased  estimator  than  the  MLE  estimator. 

Reduced  graphics  are  not  of  much  use  for  the  jacknifed 
estimators  (Figures  Sc  §  5d)  since  the  bias  is  almost  com¬ 
pletely  removed  by  jacknifing.  They  are  useful  only  in  the 
sense  that  they  are  on  the  same  scale  as  their  non jacknifed 
forms  and  a  comparison  of  bias  between  estimators  can  only 
be  done  graphically  if  the  graphics  of  estimators  are  on  the 
same  scale.  This,  however,  was  already  accomplished  from 
Figures  4a  through  4d. 

As  was  done  for  the  variance,  the  amount  of  bias  for 
a  given  subsample  size  can  be  estimated  by  dividing  the  coef¬ 
ficients,  82,  etc.,  by  the  sample  size  squared,  cubed, 
etc.,  as  appropriate  and  summing  (See  Chapter  II  "USAGE"  for 
the  mathematical  formula) .  The  computations  will  not  be 
carried  out  here  as  in  the  case  of  the  variance  but  it  is 
useful  to  compare  the  8^  coefficients  of  the  estimators 
since  these  coefficients  are  likely  to  contribute  the  most 
significant  amount  of  bias  (they  are  divided  by  the  smallest 
denominator).  In  the  case  of  the  moment  estimator,  is 
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16.50  and  for  the  MLE  estimator  it  is  14,77.  This  suspected 
greater  bias  by  the  moment  estimator  can  be  seen  visually 
from  the  graphics.  The  smaller  coefficients  for  3^  of  the 
jacknifed  moment  and  MLE  estimators  are  1.47  and  4.09  respec¬ 
tively.  This  gives  an  indication  of  how  much  less  the  bias 
becomes  when  the  estimators  are  jacknifed.  This  kind  of 
comparison  of  3^^  coefficients  should  be  done  with  caution, 
though,  since  subsequent  coefficients  (32*  @3>  ••)  niay  be  so 
large  that  they  contribute  more  to  the  bias  than  the  3^  coef¬ 
ficient.  It  should  also  be  pointed  out  that  these  3^  coeffi¬ 
cients  should  be  zero  because,  as  mentioned  earlier,  the 
jacknife  removes  the  3^(l/n)  term.  In  fact,  their  estimated 
standard  deviations  are  4,00  and  5.67  respectively,  so  that 
they  are  within  one  standard  deviation  of  zero. 

.4s  important  as  the  size  of  these  bias  coefficients 
are  their  accuracy.  Comparing  the  bias  coefficients  from 
Figures  5a  through  5d  in  relation  to  their  standard  deviation 
with  those  of  Figures  3a  through  3d  (Table  4)  which  were 
determined  from  only  two  replications,  shows  the  advantage 
of  large  replication  in  estimating  the  bias  of  an  estimator. 


51 


Table  5:  Bias  coefficients  for  each  estimator's 

population  of  estimates  and  their  standard 
deviations  for  small  and  large  replications 
(M*2  and  M*100) 


RT  AS 

NUMBER 

OF  REPLICATIONS 

ESTIMATOR 

D  X 

COEFF. 

TWO  (2) 

ONE  HUNDRED  (100) 

7.40 

16.49 

1 

(13.5) 

(2.59) 

MOMENT 

183.1 

-182.2 

(1105.5) 

(204.1) 

St 

1947.9 

3659.04 

0 

(21284.5) 

(4526.5) 

Si 

32.65 

14.77 

i 

(19.1) 

(2.83) 

MLE 

-1998.7 

-14.12 

i. 

(1031.2) 

(220.2) 

- 

St 

47308.0 

1602.74 

(18073.1) 

(4798 . 3) 

s, 

6.59 

1.47 

1 

(11.0) 

(4.01) 

JNK  MNT 

a 

1 

-912.54 

9.75 

(1016.4) 

(335.5) 

St 

21772.7 

-1369.26 

3 

(23318.6) 

(-411.7) 

Si 

75.33 

4.09 

1 

(16.23) 

(5.67) 

JNK  MLE 

^2 

-5815.6 

-203.3 

(1450.6) 

(436.02) 

112632.0 

793.11 

j 

(32037.1) 

(9243.2) 
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Figure  5d  (REDUCED  GRAPHICS) 


3 .  Case  1,  Example  5:  Large  Subsample  Sizes  and  Large 

Replications  (k»5.0) 

Figures  4a  through  4d  showed  the  necessity  of  increas¬ 
ing  subsample  sizes  in  order  to  examine  the  normality  of  the 
estimates.  Figures  6a  through  6d  are  for  increased  subsample 
sizes  (100  to  2000)  with  a  large  number  of  replications  (M»25). 
The  box  plots  in  these  figures  show  the  increase  in  normality 
with  increasing  subsample  size  n^^  when  they  are  compared  with 
small  subsample  sizes  (Figures  4a  through  4d) .  The  reduction 
in  variability  is  also  obvious  by  noticing  the  lack  of  outliers 
in  the  box  plots  and  the  smaller  range  of  the  vertical  scale 
for  larger  sample  sizes.  The  standard  devia-^ion  for  each  of 
the  estimators  for  small  subsample  sizes  (ng=500.  Figures  4a 
through  4d)  is  compared  to  those  for  large  subsample  sizes 
(ng=2000,  Figures  6a  through  6d) ,  Table  6.  It  can  be  seen 
that  the  variance  is  reduced  for  larger  subsample  sizes  by 
the  appropriate  factor  of  two.  However,  the  degrees  of  freedom 
in  the  standard  deviation  of  the  estimators  are  25  at  n^  *  2000 
and  100  at  n^  =  500.  Again  it  is  clear  from  these  figures 
that  jacknifing  causes  no  lack  of  precision  in  either  the  MLE 
or  moment  estimators. 

The  bias  in  the  moment  and  MLE  estimators  is  difficult 
to  compare  from  the  graphics  because  of  the  range  of  the  verti¬ 
cal  scale.  For  this  reason  reduced  graphics  (Figures  7a  §  7b) 
for  these  two  estimators  are  produced.  Reduced  graphics  for 
the  jacknifed  estimators  would  not  be  very  useful  since  it 


has  been  shown  in  previous  examples  that  very  little  bias 
exists.  For  this  reason  they  are  omitted. 


Table  6:  Standard  deviations  of  the  population  of 
estimates  of  the  shape  parameter  for  each 
estimator  for  subsample  sizes  500  and  2000 

SAMPLE  SIZES 


ESTIMATOR 

ng  »  500 

ng  »  20 

MOMENT 

.3417 

.1592 

MLE 

.2730 

.0910 

JACKNIFED 

MOMENT 

.3400 

.1597 

JACKNIFED 

MLE 

.  2759 

.1090 

A  change  seems  to  have  occurred  if  a  comparison  is 
made  of  Figures  7a  §  7b  with  those  of  Figures  5a  5  5b.  The 
graphics  seem  to  indicate  that  the  MLE's  unbiased  estimates 
have  shifted  from  being  biased  on  the  high  side  to  the  low 
side  for  subsample  size  nj^  =  100  and  possibly  n-,»200.  Also, 
it  appears  that  the  two  estimators  have  about  the  same  bias 
or  that  possibly  the  moment  estimator  has  less  bias.  Looking 
further  at  the  bias  coefficients  that  determine  these  graphics, 
it  can  be  seen  that  the  coefficients  in  Figures  5a  and  5b  are 
estimated  with  more  precision  than  those  of  Figures  7a  and 
7b.  One  may  be  tempted  to  assume  that  the  MLE  estimator  is 
more  biased  for  larger  sample  sizes.  But  this  is  not  likely 
and  the  change  in  the  bias  structure  for  Figure  7b  can  be 
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attributed  to  the  difference  in  the  values  of  the  bias 
coefficients  62*  ' compared  to  those  of  Figure  5b. 

These  coefficients  should  be  the  same  but  the  coefficients 
of  Figure  5b  were  computed  from  100  estimates  (100  regressions) 
and  those  from  Figure  7b  only  25.  The  precision  can  be 
expected  to  be  better  in  the  first  case.  This  is  obviously 
the  case  if  the  standard  deviation  of  the  coefficients  are 
compared,  and  more  confidence  can  be  placed  in  the  bias 
structure  of  Figure  5b. 
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B.  EXAMINING  THE  GAMMA  DISTRIBUTION  FOR  THE  SHAPE  PARAMETER 
(K=0.2) 

1 .  Case  2,  Example  1:  Small  Sub sample  Sizes  and 

Small  Replications  CK=0.2) 

In  the  preceding  case  for  K=5.0  four  estimators  of 
the  shape  parameter  of  the  gamma  distribution  were  evaluated 
for  varying  subsample  sizes  and  replications.  This  evaluation 
is  repeated  here  with  the  shape  parameter  changed  from  K=5.0 
to  K=0.2.  This  second  set  of  comparisons  is  not  only  useful 
in  demonstrating  how  the  RAGE  package  shows  the  behavior  of 
different  estimators  in  regards  to  bias  and  variance,  but 
also  hov/  these  estimators  may  differ  when  the  parameter  being 
evaluated  takes  on  different  values. 

Figures  8a  through  8d  are  graphical  results  of  the 
four  estimators  of  the  shape  parameter,  K=0.2.  For  these 
small  subsample  sizes  and  small  population  of  estimates  it 
is  apparent  from  the  graphics  that  the  MLE  estimator  (Figure 
8b)  has  tighter  box  plots  than  does  the  moment  estimator 
(Figure  8a).  This  shows  graphically  that  the  variance  for 
the  MLE  estimator  is  less  than  that  of  the  moment  estimator. 
This  was  not  obvious  when  the  shape  parameter  was  5.0.  In 
Figures  3a  §  3b  of  Case  1,  Example  1,  it  was  difficult  to 
distinguish  any  difference  in  the  box  plots.  This  comparison 
shows  how  the  MLE  estimator  is  much  better  at  reducing  the 
variance  when  the  actual  shape  parameter  is  as  small  as  0.2. 
Or,  alternatively,  it  may  be  an  indication  of  how  much  worse 
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the  moment  estimator  is  for  estimating  shape  parameters 
that  are  small. 

Figures  8c  §  8d  are  the  jacknifed  forms  of  the  esti¬ 
mators.  It  is  not  at  all  clear  from  the  graphs  that  the 
jacknife  is  affecting  the  variance  of  the  estimates.  This 
lack  of  change  could  be  attributed  to  the  small  number  of 
estimates  which  suggest  that  the  number  of  replications 
should  be  increased  to  achieve  a  better  graphical  evaluation. 

In  Case  1,  Examples  1,  2  and  3,  the  coefficients  of 
variance  (a^^,  a2 ,  aj  and  a^)  were  helpful  in  determining 
which  estimators  produced  the  least  variance.  A  look  at 
in  Figures  8b  and  8d  shows  that  they  are  negative  numbers. 

Since  often  contributes  the  most  significant  amount  to 
the  variance,  there  is  some  suspicion  that  a  determination 
of  the  variance  by  the  use  of  these  coefficients  could  result 
in  negative  variances  for  the  MLE  estimator  and  its  jacknifed 
form.  If  and  are  not  large  positive  numbers  (relatively) 
it  is  especially  likely.  In  these  two  cases  a2  is  a  rela¬ 
tively  large  positive  number  and  negative  variances  do  not 
occur  if  the  calculations  are  carried  through.  However, 
for  lesser  degrees  of  regression  (i.e.  D=2)  negative  variances 
can  occur.  This  is  a  result  of  the  wide  fluctuation  of  the 
variances  of  the  small  subsample  sizes  whose  values  are  used 
as  dependent  variables  in  the  regression  which  determine 
,  and  a-.  These  variance  coefficients  are  often 


questionable  when  the  variances  of  the  population  of  esti¬ 
mates  vary  over  a  wide  range  and  further  suggest  increasing 
the  number  of  replications  so  that  the  variance  coefficients 
can  be  more  accurately  determined. 

The  bias  in  the  moment  estimator  (Figure  8a)  appears 
to  be  considerably  more  than  that  of  the  MLE  estimator 
(Figure  8b)  when  the  graphs  of  figures  are  compared.  This 
bias  is  further  reduced  when  the  jacknife  is  applied  to  the 
estimators  (Figures  8c  8  8d) .  A  comparison  of  the  bias  is 
inadvisable  at  this  point  since  the  bias  coefficients  have 
large  standard  deviation  values  in  comparison  to  their 
estimates.  The  8^  coefficient,  for  example,  of  the  MLE 
estimator  is  .0864  and  its  standard  deviation  is  .189.  As 
mentioned  in  previous  examples,  8^  will  normally  contribute 
the  most  to  the  bias  and  since  its  standard  deviation  is 
larger  than  the  estimate  itself,  its  accuracy  is  highly 
questionable.  Before  any  firm  conclusions  can  be  drawn 
concerning  the  bias  of  the  estimators  more  estimates  are 
required,  which  implies  increasing  the  number  of  replications. 
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Case  2.  Example  2:  Small  Sub sample  Sizes  and 
Large  Replications 

Figures  9a  through  9d  are  for  small  subsample  sizes 
but  large  replications  (M=100) .  This  increased  number  of 
estimates  produces  more  definitive  graphics  for  better 
evaluations  of  the  bias  and  variance  of  the  estimates  pro¬ 
duced  by  the  estimators- 

The  graphs  of  the  moment  and  MLE  estimators  (Figures 
9a  §  9b)  show  clearly  as  before  that  the  MLE  estimator  has 
much  less  variance  than  does  the  moment  estimator.  Both, 
however,  have  positive  skewness  and  high  nonnormal  estimates 
up  to  subsample  size  n^=125.  The  jacknifed  form  of  the 
moment  estimator  has  about  the  same  variance  as  its  nonjack- 
aifed  form  but  is  more  symmetrically  distributed.  This  can 
be  seen  in  a  comparison  of  the  box  plots  of  Figures  9a  5  9c. 

The  distribution  of  the  population  of  estimates  for 
the  jacknifed  MLE  estimator  shows  little  if  any  difference 
graphically  or  statistically  than  the  nonjacknifed  estimates. 
This  can  be  seen  by  comparing  Figures  9b  5  9d ,  which  show 
little  difference  when  the  box  plots  are  compared  for 
various  subsample  sizes.  For  large  values  of  the  shape 
parameter,  i.e.  K=*5.0,  it  could  be  seen  from  Figures  4b  5  4d 
of  Case  1,  Example  1  that  the  jacknifed  MLE  estimator 
increased  the  symmetry  of  the  estimates.  This  was  obvious 
both  in  the  graphics  and  by  comparing  the  skewness  coeffi¬ 
cients.  For  K»0.2  the  only  indication  that  the  estimates 
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are  more  symmetrically  distributed  is  the  comparison  of  the 
skewness  coefficients  for  small  subsample  sizes  (n^=33, 
n2=50  and  n^*?!)  of  the  MLE  estimator  and  its  jacknifed  form 
(Figures  9b  5  9d) .  The  skewness  appears  to  decrease  slightly, 
but  not  discernably,  for  larger  subsample  sizes. 

The  amount  of  variance  that  was  observed  graphically 
can  also  be  seen  statistically  by  computing  the  variance  of 
each  estimator  for  designated  subsample  sizes  using  the  vari¬ 
ance  coefficients  and  of  each  estimator.  These 

coefficients  for  each  estimator  will  result  in  positive 
variances.  This  is  not  always  the  case  when  the  coefficients 
are  determined  from  subsample  sizes  with  a  small  population 
of  estimates.  This  e.xample  illustrates  the  increased  confi¬ 
dence  that  can  be  put  into  the  variance  coefficients  as  the 
number  of  estimates  (replications)  increases. 

The  bias  again  appears  to  be  substantially  greater 
for  the  moment  estimator  than  the  MLE  estimator  for  comparable 
subsample  sizes.  Here,  because  of  the  range  of  the  vertical 
scale,  reduced  graphics  are  produced  for  the  moment  and  MLE 
estimators  to  better  compare  the  bias  of  these  estimators 
by  reducing  the  scale  (Figures  10a  and  10b) .  A  comparison 
of  Figures  10a  and  10b  clearly  shows  the  increased  bias 
caused  by  the  moment  estimator.  It  was  shown  in  Example  1 
of  Case  1  that  the  jacknifed  estimators  have  little  bias 
remaining  in  the  estimates.  This  appears  to  be  the  case  for 
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K»0.2,  as  indicated  by  the  small  bias  coefficients  of  the 
jacknifed  estimators  (Figures  9c  §  9d) .  For  this  reason, 
reduced  graphics  for  the  jacknifed  estimators  are  omitted, 
since  indications  are  they  would  show  little  more  graphically 
than  Figures  9c  5  9d. 

It  is  interesting  to  note  that  the  bias  of  the  moment 
estimator  compared  to  that  of  the  MLE  estimator  is  not  only 
greater,  but  greater  for  small  values  of  the  shape  parameter 
(i.e.  K='0.2)  than  for  large  values  (i.e.  K=5.0).  This  can 
be  seen  by  comparing  the  bias  in  Figures  10a  §  10b  with  that 
of  Figures  5a  §  Sb.  This  indicates  that  the  moment  estimator 
is  not  only  more  biased  than  the  MLE  estimator  but  consider¬ 
ably  more  biased  when  the  shape  parameter  being  estimated  is 
small . 
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3 .  Case  2.  Example  3:  Large  Subsample  Sizes  and 

Large  Replications  (K«6723 

As  a  last  example,  graphs  for  large  subsample  sizes 
with  a  large  number  of  replications  (M*25)  are  produced  for 
each  estimator.  The  box  plots  of  these  graphs  (Figures  11a 
through  lid)  look  similar  to  their  counterparts  for  small 
sample  sizes  and  large  replications  (Figures  9a  through  9d) . 
This  is  due  to  the  difference  in  the  vertical  scale  and  does 
not  indicate  a  lack  of  variance  reduction  by  larger  sample 
sizes.  The  range  of  the  vertical  scale  of  Figures  11a 
through  lid  is  .495  and  is  smaller  than  that  of  Figures  9a 
through  9d,  which  is  .9865,  by  a  factor  or  two.  This  is  a 
simple  but  important  illustration  of  the  caution  that  must 
be  exercised  when  comparing  graphics  that  are  not  of  the 
same  scale. 

Other  than  the  reduced  variance  caused  by  larger 
sample  sizes,  the  same  conclusions  that  were  reached  in 
Example  2  of  this  case  concerning  the  variance  and  the  bias 
of  the  estimators  is  shown  to  hold  here.  The  MLE  estimator 
shows  considerably  less  variance  and  bias  than  the  moment 
estimator.  The  jacknifed  moment  estimator  results  in  more 
normally  distributed  estimates  and  less  bias.  The  jacknifed 
MLE  estimator  shows  some  improvement  toward  normality  of  its 
estimates  and  its  variance  has  changed  little. 

As  in  Case  1  for  K=5.0,  the  larger  the  subsample 
size  the  more  normal  are  the  population  of  estimates.  This 
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increased  normality  can  be  seen  graphically  by  comparing  the 
graphs  of  Figures  11a  through  lid  with  those  of  Figures  9a 
through  9d.  The  lack  of  outliers  on  the  high  side  of  the 
box  plots  for  Figures  9a  to  9d  is  an  indication  of  convergence 
normality. 

Table  7  summarizes  the  results  of  the  comparisons  of 
the  four  different  estimators  for  K=5.0  and  K=0.2.  A  com¬ 
parison  of  the  standard  deviations  shows  the  moment  estimator 
to  be  consistently  higher  than  the  MLE  estimator.  These 
standard  deviations,  however,  show  little  change  when  the 
primary  estimators  are  jacknifed,  indicating  little  or  no 
change  in  variance  by  jacknifing.  The  bias  values  of  the 
moment  estimator  are  consistently  higher  and  the  root  mean 
square  error  (r.m.s.e.)  is  considerably  larger  for  K*0.2  in 
comparison  to  the  MLE  estimator,  indicating  how  poor  the 
moment  estimator  is  for  evaluating  shape  parameters  of  gamma 
distribution  that  have  small  values. 

These  conclusions  were  easily  drawn  from  the  graphs 
and  the  statistical  output  provided  by  the  RAGE  program  and 
did  not  require  extensive  programming  and  analysis  by  the 
experimenter.  This  kind  of  expedient  examination  of  estima¬ 
tors  not  only  saves  the  experimenter  time  but  reduces  signifi¬ 
cantly  the  computer  time  required  for  analysis.  This 
facilitation  of  analysis  is  the  primary  purpose  of  the  RAGE 
methodology. 
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Table  7:  Bias  and  variance  comparisons  for  varying 
subsample  sizes  and  shape  parameters 


SAMPLE  § 

SHAPE  SUB SAMPLE 

PRMTR  SIZE  M  ESTIMATOR  ESTIMATE  STD. DEV.  BIAS  R.M.S.E. 


0.2 

500 

100 

MOMENT 

.2071 

.002927 

-.0071 

.03012 

MLE 

.2006 

.0008701 

-.006 

.01057 

JN  MOM 

.2015 

.003211 

-.0015 

.03215 

JN  MLE 

.2000 

.008707 

.00 

.  0087 

0.2 

2000 

25 

MOMENT 

.2021 

.003584 

-.0021 

.01804 

MLE 

.2001 

.0008417 

.0001 

.00421 

JN  MOM 

.2005 

. 003644 

-  .0005 

.01823 

JN  MLE 

.1997 

.0008966 

.0003 

.00494 

5.0 

500 

100 

MOMENT 

5.012 

.03417 

-.012 

.3419 

MLE 

5.095 

. 02730 

-.095 

.2891 

JN  MOM 

4.981 

.0340 

.019 

.3405 

JN  MLE 

5.057 

.02759 

-.057 

.2817 

5.0 

2000 

25 

MOMENT 

4.994 

. 03184 

.006 

.1593 

MLE 

5.049 

.01821 

-  .049 

.1033 

JN  MOM 

4.988 

.03194 

.  012 

.16015 

JN  MLE 

5.014 

.02180 

-  .014 

.1099 
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